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ABSTRACT 

Protoplanetary disks are believed to accrete onto their central T Tauri star because of magnetic stresses. 
Recently published shearing box simulations indicate that Ohmic resistivity, ambipolar diffusion and the Hall 
effect all play important roles in disk evolution. In the presence of a vertical magnetic held, the disk remains 
laminar between 1-5 au, and a magnetocentrifugal disk wind forms that provides an important mechanism for 
removing angular momentum. Questions remain, however, about the establishment of a true physical wind 
solution in the shearing box simulations because of the symmetries inherent in the local approximation. We 
present global MHD simulations of protoplanetary disks that include Ohmic resistivity and ambipolar diffusion, 
where the time-dependent gas-phase electron and ion fractions are computed under EUV and X-ray ionization 
with a simplihed recombination chemistry. Our results show that the disk remains laminar, and that a physical 
wind solution arises naturally in global disk models. The wind is sufficiently efficient to explain the observed 
accretion rates. Eurthermore, the ionization fraction at intermediate disk heights is large enough for magneto- 
rotational channel modes to grow and subsequently develop into belts of horizontal held. Depending on the 
ionization fraction, these can remain quasi-global, or break-up into discrete islands of coherent held polarity. 
The disk models we present here show a dramatic departure from our earlier models including Ohmic resistivity 
only. It will be important to examine how the Hall effect modihes the evolution, and to explore the inhuence 
this has on the observational appearance of such systems, and on planet formation and migration. 

Subject headings: accretion, accretion disks - MHD - methods: numerical - protoplanetary disks 


1. INTRODUCTION 

Understanding the complex dynamical evolution of proto¬ 
planetary disks (PPDs) is of key interest both for building a 
comprehensive theory of planet formation, and as a means of 
explainin g the observationall y inferred properties of these ob¬ 
jects (see Turner et al.|2014| for a recent review). Eor exam¬ 
ple, PPDs are known to accrete gas onto their host stys at a 


typical rate of Mq yr“^ (|Gullbring et al. 1998 Hart 

|mann et al.|19981l and have li fe tirnes in the range 3-10 Myr 
( Sicilia-Aguilar et al.||2006 i. More recent observations have 
indicated the presence of turbulence in the disks surrounding 
the young stars TW Hydrae and H D 163296, based on a naly- 
sis of their molecular line emission (|Hughes et al.|2011||. Eur- 
ther evidence for turbulent broadening comes, for insta nce, 
from the CO ro vibrational band head shape measured by |Na-| 
[jita et al.|d2009| l in V1331 Cygni. During the T Tauri (class II) 
phase, self-gravity is expected to provide a negligible contri¬ 
bution to angular momentum transport (due to the low disk- 
to-star mass ratio), and disk accretion is instead believed to 
be driven by magnetic effects. Two possible sources of angu¬ 
lar momentu m transport are through a magnetocentrifugally 
driven wind ( |Blandford & Payiie] 1982[ [Pudritz & Norman] 
|1983|l, or through the operation of the rnagnetorotational in¬ 
stability (MRI,|Balbus & Hawley|1991|l, whose nonlinear out¬ 
come in the ideal-MHD limit is niagnetohydrodynamic turbu¬ 


lence ( [Hawley & Balbu^|1991[ [Hawley, Gammie, & Balbils 

Except for the innermost regions of PPDs, where the tem¬ 
perature T 1000 K, and the ionization of alkali metals such 

oliver.gressel@nbi.dk (OG); neal.turner@jpl.nasa.gov (NJT); 
r.p.nelson@qmul.ac.uk (RPN); cmcnally@nbi.dk (CPMN) 


as sodium and potassium pr ovides strong coupling between 
the gas and magnetic held (jUmebayashi & Nakano||1988|), 
non-ideal MHD effects such as Ohmic resistivity, arnDipoto 
diffusion (AD) and Hall drift are expected to be important due 
to the low ionization levels of the gas (e.g. [Blaes & Balbus 


1994 


Wardle|1999[[Sano, Miyama, Umebayashi, & Nakano 

Balbus & Terc[uem|2001|l. External sources of ioniza- 


tion such as stellar X-rays, EU V photons and galactic cosmic 
rays are expected to ionize the disk surface layers, provid¬ 
ing good coupling there - even though recent results suggest 
that cosmic rays may be highly attenuated by the stellar wind 
( jCleeves et al.|2013] l. Nevertheless, deep within the disk the 
evolution should be co ntrolled by the n on-ideal effects. The 
recognition of this led [Gammiej (|1996|l to propose what has 
now become the traditional dead-zone model in which the 
disk surface layers accrete by sustaining MRI turbulence, with 
the shielded interior maintaining an inert and magnetically de¬ 
coupled “dead zone”, where the MRI is quenched by the ac¬ 
tion of Ohmic resistivity only. The potential importan ce of the 
other non-ideal effects has long been recognized (e.g.jSano &| 
Stone|2002| jSalmeron & Wardle|2003|l, but it is only recently 
that nonlinear shearing box simulations including ambipolar 
diffusion and the Hall term have been performed in relevant 
parameter regimes, leading to a modified picture of how disks 
accrete that deviates signihcantly from the traditional dead 
zone one. 

Generally speaking, it is expected that in disk regions be¬ 
tween 1-5 au. Ohmic resistivity will be dominant near the 
midplane, the Hall effect will be dominant at intermediate 
disk heights and ambipolar diffusion will be most important 


in low density regions higher up in the disk (Salmeron & War- 
dle|2008 1. The highest altitude surface layers are expected to 
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be ionized by stellar FUV photons, and as such will e volve in 
the ideal MHD limit ( Perez-Becker & Chiang|20TT]l . Shear¬ 
ing box simulations presented by|Bai & Stone|(|2013|l that in¬ 
clude Ohmic resistivity and ambipolar diffusion for models 
computed at 1 au demonstrate that AD quenches MRI tur¬ 
bulence. They conclude that disks will be completely lam¬ 
inar between 1-5 au with angular momentum transport oc¬ 
curring because a magnetocentrifugal wind is launched from 
near the disk surface. For reasonable assumptions about the 
magnetic field strength and geometry, accretion rates in agree¬ 
ment with the observations are obtained. The inclusion of 
the Hall term in the presence of a vertical magnetic field in¬ 
troduces dichotomous behavior, arising from the fact that the 
coupled dynamics depends on the sign of 17 • B, where 17 and 
B signify the disk angular momentum vector and the ambient 
magnetic field dire ction, respectively. She aring box s imula- 
tions presented by [Lesur et'af] ( |2014| l and |Bai ( |2014| l show 
that when 17 • B > 0, the Hall effect leads to the amplifi¬ 
cation of the horizontal magnetic fields within approximately 
two scale heights of the disk midplane, and the generation of 
a large scale Maxwell stress through magnetic braking. In ad¬ 
dition, the magnetocentrifugally driven wind is also seen to 
strengthen. When 17 • B < 0, however, magnetic stresses 
and winds are seen to weaken relative to the opposite case, 
and relative to the evolution observed when Hall effects are 
neglected. Quite how this Hall dichotomy maps onto observa¬ 
tions of PPDs remains an open question. 

The internal dynamics of protoplanetary disks in the re¬ 
gion 1-5 au are clearly of great significance for planet for¬ 
mation. The growth and settling of grains depends on the 
level of turbulence, with a laminar disk or one displaying 
weak turbulence providing the most favorable conditions - 
although it should be noted that some turbulent mixing is re¬ 
quired to maintain a population of small grains in the sur¬ 
face regions of PPDs so that their spectral energy distribution 
can be reproduced by radiative transfer models ( |Dullemond| 
& Dominik||2005 1. The collisional growth of planetesimals 


has be en'sEown to be affected str ongly by the level of turbu¬ 
lence (Gressel et al. |2011 2012[ ), and the migration of low 
mass planets is also highly sensitive to the presence or oth¬ 
erwise of turbulence, with significant st ress levels being re¬ 
quired to unsaturate corota tion torques ( |Paardekooper et ah] 
|201 1 [Baruteau et al.|20iT| l. 

This paper is the latest in a series in which we are exploring 
the influence of magnetic fields and non ideal MHD effects 
on the formation of planets, with the eventual goal of produc¬ 
ing comprehensive models of PPDs that will be us ed to study 
planet building and evolution. In ear lier work ( [Nelson &| 
|Gressel|2010[[Gressel et al.|201 ij[2012[ ) we have used a com¬ 
bination of global and shearing box simulations to examine 
the dynamical evolution of dust grains, boulders and planetes¬ 
imals in turbulent disks, with and without dead zones. More 
recently, we have studied the influence of magnetic fields on 
gap formation and gas accretion onto a forming giant planet 
using global simulat ions that also included Ohmic resistivity 
( Gressel et al.[[2013] l. In this paper, we include both Ohmic 
resistivity and ambipolar diffusion in global disk simulations, 
and follow the dynamical evolution of the resulting disk mod¬ 
els as a precursor to examining how ambipolar diffusion af¬ 
fects gas accretion onto a giant planet. Of particular interest 
is the nature of the accretion flow in the disk, the nature of any 
magnetocentrifugal wind that is launched, and how these vary 
as a function of small changes in model parameters. These 
are the first quasi-global simulations to be conducted of PPDs 


that are threaded by vertical magnetic fields and which in¬ 
clude this combination of non ideal MHD effects, and as such 
are most comparable with the shearing box sim ulations pre¬ 
sented by Bai & Stone! (20131, hereafter BS13 A fundamen¬ 
tal question raised by the shewing box simulations is whether 
or not a physical wind solution with the correct held geom¬ 
etry emerges spontaneously in global disk simulations in a 
way that is not generally observed in shearing boxes because 
of their radial symmetry. Probably the most important result 
contained in this paper is that physical wind solutions do in¬ 
deed arise spontaneously in our global simulations, demon¬ 
strating the robustness of m any of the features obtained in the 
shearing box simulations of[BS13[ 

This paper is organized as follows: In Section 1^ we de¬ 
scribe the numerical method, the disk model, as wml as the 
ionization chemistry. Our results are presented in Section!^ 
where we mainly focus on the emerging wind solutions and 
the dynamical st abili ty and evolution of forming current lay¬ 
ers. In Section [3.7| we will moreover report an instability 
that is driven by the sharp transition into the FUV dominated 
layer, and in Section [US we assess whether secondary insta¬ 
bilities can drive non-axisymmetric evolution. We summarize 
our results in Section and discuss implications for planet 
formation theory in Section]^ 


2. METHODS 

We present magnetohydrodynamic (MHD) simulations of 
protoplanetary accretion disks employing 2D axisymmetric 
and 3D uniformly-spaced spherical-polar meshes. In the fol¬ 
lowing, the coordinates (r, 6, (p) denote spherical radius, co¬ 
latitude and azimuth, respectively. We moreover ignore ex¬ 
plicit factors of the magnetic permeability po in our nota¬ 
tion. Simulations were run using the single-fluid NIRVANA- 
III code, which implem ents a second-order-accurate Godunov 
scheme ( Ziegler[[2004[ l. The code is f ormulated on orthogo¬ 
nal curvilinear rneshes ( [Ziegler[[2011|l and emplo ys the con¬ 
strained transport method ^Evans & Hawley[1988[ l to maintain 
the divergence-free constraint of the magnetic held. As an al¬ 
teration to the publicly available distribution of t he code, we 
here adopt the upwind reconstruction technique of Gardiner & 


Stone ( [2008| l to obtain the edge-centered electric field compo¬ 
nents needed for the magnetic field update (also see [Skinner[ 
[& Ostriker[2010[ l. This modification became necessary to take 
advantag e of the more accurate HLL D approximate Riemann 
solver of[Miyoshi & Kusano[([2005[l, which offers advanced 
numerical accuracy for flows in the subsonic regime. 

Our numerical model is based on solving the standard re¬ 
sistive MHD equations including Ohmic resistivity as well as 
an electromotive force resulting from the mutual collision of 
ions and neutrals. Given the typical number densities within 
PPDs, the applicability of the strong-couplin g approxim ation 
is warranted by detailed chemical modeling ( Bai|201 1} . Ac¬ 
cordingly the gas dynamics can be modeled by a single-fluid 
representing the motion of the neutral species, and the ad- 
ditio nal term can simp ly be evolved in a time-explicit fash¬ 
ion ([Choi et al.[[2009|l. The total-energy formulation of the 


NIRVANA-III code is expressed in conserved variables p, pv, 
and e = e -|- pv^/2 + B^/2, and if we define the total pressure, 
p*, as the sum of the gas and magnetic pressures, we can write 
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our equation system as 

atp +V-(pv) = 0, 

dt (pv) + V • [pvv — BB] = —pV $ , 

dte + V- [(e + p*)v - (vB)B] = V • [ (Eq +Ead) x B] 

-p(v<i>)-v + r, 

9tB — V X (V X B + Eq + Ead ) = 0) (1) 

where the electromotive forces stemming from the Ohmic and 
ambipolar diffusion terms are given by 


Eo = —Vo (VxB). 


and 


EaD = PAD 


(VxB) X b X b , 


( 2 ) 

(3) 


(with b = B/|B| the unit vector along B), respectively. The 
gravitational potential <i)(r) = —GAIq/v is a simple point- 
mass potential of a solar-mass star. We moreover obtain the 
gas pressure as p = (7 — l)e, where 7 = 7/5 is appropriate 
for a diatomic molecular gas. The source term T, mimicking 
radiative heating and cooling, is included to relax the specific 
thermal energy density e to its initial radial temperature pro¬ 
file. The relaxation is done on a fraction of the local dynam¬ 
ical time scale, which is short enough to avoid strong devia¬ 
tions from the equilibrium model but long enough to suppress 
vertical corrugation o f the disk due to the Goldreich-Schub ert- 
Fricke instability (see [Nelson, Gressel, & Umurhan|2013| for 
a detailed study of this phenomenon). We remark that, while 
this instability is physical in nature, its appearance may be ex¬ 
aggerated in a strictly isothermal simulation, and more realis¬ 
tic modeling including radiative transfer should be employed 
to study its relevance. 

Because we use a time-explicit method, large peak val¬ 
ues in the dissipation coefficients rjo and pad impose severe 
constraints on the numerically permissible integration time 
step. We address this by using a state-of-the-art s econd order 
super- time-stepping scheme recently proposed by |Meyer etal.| 
([ 20 T 2 ]|. In comparison with conventional first-order methods, 
this Runge-Kutta-Legendre scheme is free of adjustable pa¬ 
rameters and has been proven superior in terms of both accu¬ 
racy and robustness. To maintain the second-order accuracy 
of our time integration, we use Strang-splitting for the diffu¬ 
sive terms. 


In our previo us local simulations (cf. appendix B1 in |Gres- 
sel et al. |2012 1 , we have found that applying a constant cap 
on po can qualitatively change the way in which the top and 
bottom disk layers are coupled to each other. For a typical 
wind configuration, the horizontal field components generally 
change sign at the midplane. We imagine that in this situa¬ 
tion a clipped po would affect the amount of field diffusing 
into the midplane region and conversely the strength of the 
current sheets above and below the magnetically decoupled 
layer. We therefore choose not to apply a cap on po. How- 
eve r, we li mit pad such that Aad /3p = 2 Rm ad >0.1 (also 
see BSD 1 , where /3p = 2p/B^ is the plasma parameter, and 


Rm ad = Cs iT/PAD is the equivalent of a magnetic Reynolds 
number for AD. This we find greatly reduces the numerical 
cost without noticeably altering the obtained solution. 

2.1. External ionization and disk chemistry 

The Ohmic and ambipolar diffusivities are evaluated dy¬ 
namically for each grid-cell and are based on a precomputed 


look-up table, that has been derived self-consistently from a 
detailed chemical model accounting for grain-charging and 
gas-phase chemistry. The ionization rates that enter this equi¬ 
librium chemistry are governed by ionizing radiation entering 
the disk. The column densities that govern the attenuation of 
the external ionizing radiation are also computed on-the-fly, 
that is, they reflect the instantaneous state of the disc. In prin¬ 
ciple, this allows for a self-lim iting of the emerging w ind via 
shielding of ionizing radiation (|Bans & K6nigl|20121l. 

In addition to some minor modifications to the grain charg¬ 
ing prescription (as detailed below), we have improved the re¬ 
alism of our ionization model compared to our previous work. 
The main alterations concern the inclusion of two additional 
radiation sources, an un-scattered direct X-ray component, 
and hard UV photons. These have in common a shallow pen¬ 
etration depth but comparatively high flux levels, thus mainly 
affecting the ionization level of the disk’s surface layers. 


2.1.1. FUV ionization layer 


Adopting a very simple prescri 

ption based on the recent 

study by Perez-Becker & Chiang ( 

201111 , we have added the 

effect of an FUV ionization layer 

based on an assumed ion- 


ization fraction of / = 2 x 10 “® (representing full ionization 
of the gas-phase carbon and sulfur atoms susceptible to losing 
electrons when struck by FUV photons), and a collision rate 
coefficient of 2 x 10“® s“^. We further assume that FUV 

photons penetrate to a gas column density of 0.03 g cm“^. 
Due to their lower amplitude, we ignore any scattered, dif¬ 
fuse or ambient FUV illumination of the disk and only evalu¬ 
ate sight-lines pointing directly towards the cen tral sta r. Note 
that this deviates from the local simulations of |BS13| where 
the vertical gas column was used for attenuation of the inci¬ 
dent flux. Our treatment is motivated by the assumption that 
a large fraction of the FUV photons originate dir ectly from 
the central star (see upper right panel of figure 9 in|Bethell & 
Bergin|201 1[ who study the forward scattering of FU V pho¬ 
tons in detail). 


2.1.2. Illumination by X-rays 

A similar treatment is adopted for the unscattered X-ray 
component, for which w e use the fit to the Mo nte-Carlo 
radiative-t ransfer results of Igea & Glass gold ( 1999} given by 
eq. (21) in|Bai & Goodman] (|2009jl. Deviating froiii our pre¬ 
vious local simulations (cf eg (1) in |Gressel et al.|2011| l, our 
new X-ray illumination is 


CxR = 10 S ^ 


exp - 


-r^r.exp-rM 


Esc / 


Esc / 


- 6 x 10 


- 12 g-l 


exp — 


Eq 

Eab 


„-2 


(4) 


where Ea (Eb) is the gas column to the top (bottom) of 
the disk surface, and Yjq is the column density along radial 
rays towards the star. The radial column contains a contri¬ 
bution from the inner disk (which is not part of the com¬ 
putational domain) and reaches in to an in ner truncation ra¬ 
dius o f rtr = 5 i? 0 . Adopting values from|Bai & Goodman 
(20091, the shape coefficients are Esc = 7 x 10 ^^ cm“^, and 
a = 0.65 for scattered X-rays, and Eab = 1-5 x 10^^ cm“^, 
and /? = 0.4, for absorption of the direct component, respec¬ 
tively. Both contributions additionally decay with the squared 
radius to account for the decrease in X-ray luminosity. To 
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Elsosser numbers at r=1.0ou 



height z [H] 


Figure 1. Profiles of the initial Elsasser numbers (at radius r = 1 au) for 
ohmic resistivity (light orange), and ambipolai* diffusion (red). Solid lines 
represent the fiducial case of a dust-to-gas mass ratio of 10“^, whereas dash- 
dotted lines show profiles for a value of 10“^. We also plot the initial profile 
of the plasma (3 parameter (solid black line). 


account for typical median values of 10 ^° erg in observa¬ 
tional data of lum inosities in young sta r clusters such as in 
the Orion Nebula ( jGarmire et al.|[2000| l, we enhance the in¬ 
cident X-ray flux by a factor of 5x compared to the above 
stated values. In view of future work, we envisage additional 
improvements to the X-ray ionization model by adopting the 
new results of |Ercolano & Glassgold| ( |2013| l, who consider 
the enhanced ionizing effects of X-rays due to the presence 
of heavier elements (assuming solar abundance). These alter¬ 
ations have been found to lead to enhanced ionizing fluxes at 
in termediate column densiti es compared to the original results 
of |Igea & GlassgoId| ( |1999| ). 

2.1.3. Modifications to the disk chemistry 

Our chemic al network is based on the one used by lllgner"^ 
|Nelson| ( [2006 1 , labeled mode 14, with grain surface reactions 
and a simplified gas-phase reaction set involving one repre¬ 
sentative metal and one molecular ion. Small changes in re¬ 
cent years include correc ting the el ectron sticking probabili¬ 
ties for the grain charge ( |Bai|201 f] ) and consistently treating 
the molecular ion such as HCO"^ (with a mass of 29 protons) 
since this is the long-lived species in a series of several reac¬ 
tions. Further details on the re action network and d iffusivities 
can be found in section 2.2. o f Landry et al. ( 2013[ l as well as 
section 4.2 of|Mohanty et al.|(|2013|l. 

The resulting inif/a/ ionization profile expressed in Elsasser 
numbers Aq/ad = i^Vo/AD)~^ with vaz = is 

plotted in Fig. at a radius of 1 au, together with the height- 
dependent plasma /3 parameter. The impl icati ons of this figure 
for our simulations are discussed in Sect |3.1| The region with 
\z\^5 H, where the plasma parameter becomes constant with 
height is caused by applying a lower limit on the gas density 
to avo id excessive Alfven speeds in the disk halo. Following 
BSD we chose this limit at 10“® (in our global model, this 


value is with respect to the initial midplane value at any given 
cylindrical radius). We remark that the equilibrium density 
profile is artificially steep in our isothermal disk model with 
temperature constant on cylinders. In this respect, the density 


floor acts to mimic a disk with realistic radiation thermody¬ 
namics, where the irradiation heating of the upper disk layers 
naturally creates an increased pressure scale height and hence 
a much shallower density profile. Because pad (and hence 
Aad) depends on the magnetization of the plasma, the actual 
density profile in the disk halo has an influence on how well 
the gas is coupled to the magnetic held lines, which in turn af¬ 
fects the efficiency of the magneto-centrifugal wind launching 
mechanism. This caveat should be kept in mind when inter¬ 
preting the mass loading of the wind, which we expect to be a 
function of the thermodynamics in real protoplanetary disks. 


2 .2. Global disk model 

We now describe the underlying protostell ar disk model, 
which is largely identical to the one used in |Gressel et al. 
( 2013[ ). The equilibrium disk model is based on a locally- 
isothermal temperature, T, decreasing with cylindrical ra¬ 
dius as T{R) = To {R/RqY. For q = —1, such a profile 
leads to a constant opening angle throughout the disk as it 
is commonly used for global numerical simulations of PPDs, 
and which provides us with the reference disk model here. 
To study what effect disk flaring has on the absorption of 
direct FUV photons, we also produce a moderately flaring 
disk with q — —3/4. For this value, the resulting power- 
law index for the local opening angle, h{R) = H{R)/R, is 
Ip = {q + l)/2 = 0.125, which is in agreement with observa- 
tional constraints for th is value based on sub-mm observations 
( Andrews et al.||2009| l, and from multi-wavelength imaging 
( Pinte et al.|2008 i ^ The radial power-law exponents for the 
gas surface density, E(i?), are — 1/2 for our fiducial model, 
and —3/8 for the flaring disk model, respectively. 

Our temperature distribution is complemented by a power- 
law for the midplane density, PmidiR) — Po {R/RoY^ with 
p = —3/2. The equilibrium disk solution is given by 


, . . RV ( gm^ 

Xr)=«(5j-1 “P( — 


1 1 

r R 


0 (r) = Ok(T) 


(p + 9) ( f j + (1 + g) - ^ 


and (5) 

, (6) 


which can be derived from the requirement of hydrostatic/- 
dynamic force balance in the vertical and radial direc¬ 
tions, and where the Keplerian angular velocity ilyi{R) = 
yjG Mq R~^/'^. In the above equations, the square of the 
isothermal sound speed results from the prescribed tempera¬ 
ture profile as Cg = c^q {R/RqY, with a parameter Cso corre¬ 
sponding to a value of H{R) = = 0.05i?. 

Guided by previous results of local 3D simulations that ex¬ 
hibit quasi-axisymmetric structure, we mainly focus on 2D- 
axisymmetric simulations. Our computational domain covers 
r G [0.5, 5.5] au. In the latitudinal direction, the grid covers 
eight pressure scale heights on each side of the disk midplane, 
that is, 0 € [ 7 r /2 — '( 9 , 7 r /2 + i?], with d — 8H/r. In the 
case of the flaring disk, the coverage of scale heights varies 
from 8H at r = 1 au to 6.5 H at r = 5 au. The grid reso¬ 
lution is Nr X Ng X Nff, = 512 X 192 x 64, which means 
that we are able to resolve relevant MRI wave lengths. In the 
axisymmetric models we use Nr x Ng = 1024 x 384 cells, 
corresponding to > 24 grid points per pressure scale height. 


* At least when assuming well-mixed dust grains, since both observations 
are sensing the dust content of the disk rather than the gas density. 
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matching the res olution of the three-dimensional box simu¬ 
lations of|BS13| For the two non-axisymmetric simulations, 
the azimuthal extent is restricted to (j> G [0, ip], with ip = 7r/2 
(a quarter wedge), or 7r/4 to limit the computational expense. 
We furthermore note that in the ideal-MHD case, where tur¬ 
bulence develops, the azimuthal domain size has been found 
to have an effect on the final saturated state ( |Beckwith et al.| 
|201 1[ [Flock et al.|2012[ jSorathia et al.|201^ , at least for the 

case of zero net flux. 

Because our ionization model depends on real physical val¬ 
ues of the gas column density, we have to chose a meaningful 
normalization factor p^, which we fix according to a surface 
density of E = 150gcm“^ at the location R = 5au, com- 
patible with the typical minimum-mass solar nebula (MMNS 
Hayashi||l981|l. In physical units, the disk temperature is 


T = 540 K ait a radius of i? = 1 au, and T = 108 K at 
5 au, resembling typical expected conditions within the pro¬ 
tosolar disk. While our values for p, T and E are similar to the 
MMSN values at 5 au, our adoption of different values for the 
radial power-law indices means that t hese va lues are different 
at 1 au compared to those adopted by|BS13| 

Our model disk is initially threaded by a weak vertical 
magnetic field B^^R) corresponding to a midplane /3po = 
‘2p/B^ = 10^ independent of radius for our fiducial disk. 
To achieve this, the vertical magnetic flux has to vary as a 
power-law in radius, taking into account the radial distribu¬ 
tion of the midplane gas pressure that itself depends on the 
density and temperature. In our disk model, the gas pressure, 
p{R) decreases outward, implying that Bz{R) falls off with 
radius, too. While such a configuration is generally expected 
when accounting for inward a dvection and outward diffusion 
of the embedded vertical flux ( Okuzumi et al.|20T4| , our par¬ 
ticular choice of keeping the relative field strength constant 
is of course arbitrary. To preserve the solenoidal constraint 
to machine accuracy, we compute the poloidal field compo¬ 
nents from a suitably defined vector potential A^{r,9). The 
initial disk model is perturbed with random white-noise fluc¬ 
tuations in the three velocity components. The rms amplitude 
of the perturbation is chosen at one percent of the local sound 
speed. We furthermore add a white-noise perturbation in the 
magnetic field with an rms amplitude of a few pG, which is 
on the sub-percent level compared with the net-vertical field 
of typically a few ten mG. 

2.3. Boundary conditions 

To complete the description of our numerical setup, we 
need to specify boundary conditions (BCs). For the fluid 
variables, we employ standard “outflow” conditions at the in¬ 
ner and outer radial domain edges, rjn and rout- This type 
is a combination of a zero-gradient condition in the case of 
n • v(rin, 9) > 0 (where n denotes the outward-pointing nor¬ 
mal vector), and reflective BCs in the opposite case, thus pre¬ 
venting inflow of material from outside the domain. At the 
upper and lower boundaries (that is, in the 9 direction), we fur¬ 
thermore balance the ghost zone values such that a hydrostatic 
equilibrium is obtained. This helps to control artificial jumps 
of the fluid variables adjacent to the domain boundary as they 
are frequently encountered with finite-volume schemes. 

In contrast to our earlier global simulations of layered ac¬ 
cretion disks subject to Oh mic resistivity and con taining an 
embedded gas-giant planet (jGressel et al.||2013|l, we here 
make a different choice for the vertical magnkic-field bound¬ 
ary condition. Whereas we previously applied magnetic BCs 
of the perfect conductor type (that is, forcing the normal field 


component to zero and allowing non-zero parallel field), we 
here use pseudo vacuum conditions, conversely enforcing the 
perturbed parallel field to vanish at the surface and only allow¬ 
ing a perpendicular perturbed field. Note that we exclude the 
initial net-vertical field from the procedure such that only de¬ 
viations from the background field are subject to the normal- 
field condition. The vertical flux threading the disk is thus 
preserved. Since the disk’s upper layers are magnetically 
dominated, and the Lorentz force acts to align the flow lines 
with the magnetic field, letting the field lines cross the domain 
boundary is of course a prerequisite for launching an outflow. 
While we realize that forcing the radial and azimuthal com¬ 
ponents of the perturbed field to vanish at the boundary may 
unduly restrict the magnetic topology of the emerging wind 
solution, we postpone the study of less-restrictive but more 
cumbersome boundaries to future work. 

Unlike in a radially-periodic local shearing-box simulation, 
our global model is critically affected by the inner radial 
boundary condition. In a real protoplanetary disk, we can 
expect that alkali metals will be thermally ionized in this re¬ 
gion, leading to the development of the MRI on timescales 
short compared with the orbital timescale at the inner domain 
boundary of our model. This poses the question how to best 
interface the MRI-turbulent inner disk with the magnetically 
decoupled midplane of the outer disk that we model. It is 
likely that MRI-generated fields can efficiently leak into the 
outer disk via magnetic diffusion. In a first attempt to account 
for the MRI-active inner disk, we gradually taper the diffu- 
sivity coefficients to zero within the inner 0.5 au of our disk 
model. Studying the inner edge of the dead zone will however 
re quire dedicated simulatio ns (similar to the ones performed 
by Dzyurkevich et al.|2010) l including this transition region. 

3. RESULTS 


The main ai m of o ur paper is to establish the laminar wind 
solutions that |BS13| previously found in local shearing box 
simulations in the context of global disk simulations. In the 
interest of economic use of computational resources and to 
guarantee adequate numerical resolution of our global mod¬ 
els, we primarily focus on 2D axisymmetric simulations, but 
have also run three-dimensional simulations to check for non- 
axisymmetric solutions. Since all our runs include a net- 
vertical flux, axisymmetry is warranted to obtain a reading 
on the development of the MRI. In cases where the solution 
proves laminar, axisymmetry is likely to produce a reason¬ 
ably accurate picture. We address the possibility of non- 
axisymmetric secondary instabilities using a limit ed s et of 
three-dimensional calculations, described in Section jUSl 
The simulations conducted for this work are listed in Ta¬ 
ble where we summarize key model parameters. We as¬ 
sume a fiducial dust-to-gas mass ratio of 10“^, that is, a de¬ 
pletion by a factor of ten compared to interstellar abundances. 
Starting from the classic layered PPD (model ‘0-b6’) includ¬ 
ing only Ohmic diffusivity, we compare this standard “dead- 
zone” disk with a corresponding model, ‘OA-b6’, additionally 
including the effects of ambipolar diffusion. Our fiducial ref¬ 
erence model is ‘OA-b5’, with a midplane plasma parameter 
of 10^. In the presence of combined ambipolar diffusion and 
Ohmic resistivity, we observe the waning of the MRI, which 
is instead replaced by a laminar wind solution. Unlikem ge¬ 
ometrically constrained shearing box simulations (|BS13| [Bat] 
2013| l, we naturally obtain a field configuration with field lines 


bending outward on both “hemispheres” of the disk. Notably, 
this topology produces thin current layers, which have previ- 
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Table 1 

Summary of simulation runs. 


Run label 

Ohm 

AD 

/3p0 

d/g 

<? 

Resol. 

0-b6 

o 

- 

“TcF" 

10“® 

-1 

1024x384 

OA-b5 

o 

o 

10® 

10“® 

-1 

1024x384 

OA-b6 

o 

o 

10® 

10-® 

-1 

1024x384 

OA-b7 

o 

o 

107 

10-® 

-1 

1024x384 

OA-b5-d4 

o 

o 

10® 

10-4 

-1 

1024x384 

OA-b5-flr 

o 

o 

10® 

10“® 

-3/4 

1024x384 

OA-b5-flr-nx 

o 

o 

10® 

10“® 

-3/4 

512x192x128 

OA-b5-nx 

o 

o 

10® 

10“® 

-1 

512x192x64 


Models are labeled according to the included micro-physical effects (first two 
letters) and the strength of the net-vertical magnetic field (expressed in terms 
of the midplane value /3p o, prefixed with the letter ‘b’). Further labels refer to 
the dust-to-gas mass ratio (‘d/g’, prefixed with the letter ‘d’), the power-law 
index, q, of the radial temperature profile (‘fir’ for “flaring”), and whether the 
azimuthal dimension is included (‘nx’ for “non-axisymmetric”). 


ously been discussed by |BS13[ The stability and evolution of 
these current layers will be one focus of our paper. 

We perform additional analysis on the influence of further 
key input parameters. With the global disk model allowing 
direct illumination from the central star, it becomes possible 
to address the question of how the disk’s ionization is affected 
by disk flaring. This is studied in model ‘OA-b5-flr’, where 
we use a power-law index of q = —3/4 instead of q = —1 
to obtain a m oderately flaring disk surface as is suppor ted by 
observations ( jPinte et'aflpOOSl [Andrews et al.||2009| l. The 
role of dust depletion, driven by processes such as coagula¬ 
tion into larger grains and/or sedimentation, is considered in 
model ‘OA-b5-d4’ with a reduced dust-to-gas mass fraction 
of 10“'* (as used in BS13| compared with our fiducial value 
of 10“^. For the sake of brevity, we here refrain from varying 
any of the many other input parameters of our model, as for 
instance, the X-ray or CR intensities, or the FUV penetration 
depth, as well as parameters affecting the disk thermodynam¬ 
ics. 


3.1. Preliminaries 

The in puts to our models are very similar to those included 
in |BS13| so it is instructive to compare the Elsasser number 
profiles in our model with their fiducial model as a means of 
understanding the similarities and differences betw een the ir 
results and ours. The fiducial model presented in |BS13| is 
computed at 1 au, and has a dust-to-gas mass fraction of 10 *, 
so we can compare this with the dot-dashe d lines i n Figure[T] 
The profiles shown there, and in figure 1 of |BS13| are similar, 
with /3p decreasing below unity at disk heights \z\ > 4.57T, 
preventing MRI turbulence from developing at these high alti¬ 
tudes. In their initial condition and ours, Aq increases mono- 
tonically from the midplane upward. The two initial states 
also have similar profiles for Aad. displaying local maxima 
at intermediate disk heights ( z ^ zL 2.5H). The differences 
in the surface densities in the |BS 13] models and ours at 1 au, 
combined with our inclusion of a direct X-ray component in 
addition to the scattered component means that Aad 10 at 
this local maximu m in ou r model, whereas it only rises mod¬ 
estly above unity in BS13 We therefore expect that our model 
with a dust-to-gas mass fraction of 10“^ should contain nar¬ 
row regions at intermediate heights that s upport the growth of 
MRI channel modes. It is noteworthy that|BS13|also observed 
the development of the MRI at the beginning of their simula¬ 
tions, but report that the subsequent amplification of the held 


causes the ambipolar diffusion to increase. Their disk then re¬ 
laxes to completely laminar state in which angular momentum 
transport is driven by a magneto-centrifugal wind. The larger 
value of Aad in our models may allow MRI turbulence to be 
sustained in these regions, or may instead lead to quenching 
of the MRI after it has reached a more nonlinear stage of de¬ 
velopment. In this paper, we define our fiducial model to be 
one in which the dust-to-gas mass fraction is 10“^, and the 
Elsasser numbers for this case are shown by the solid lines in 
Fig. 0 The larger dust concentration reduces the gas-phase 
electron fraction and Elsasser numbers, and the local max¬ 
imum in Aad now peaks moderately above unity ( but stil l 
attaining a larger peak value than the fiducial model in BS 131. 


3.2. Traditional dead zone model 

We begin discussion of the simulation results by consider¬ 
ing the 2D axisymmetric run 0-b6, for which Ohmic resis¬ 
tivity is included and ambipolar diffusion is neglected. To 
introduce our setup, and give the reader an impression of the 
general appearance of our PPD model. Pig. S visualizes the 
magnetic field and poloidal flow vectors from me model. The 
color coding of the azimuthal field strength shows the MRI- 
turbulent surface layers with tangled poloidal field lines super¬ 
imposed in white. In the outer disk, where the MRI has not 
fully developed yet, the parity of the horizontal field is odd 
with respect to reflection at the z = 0 line, consistent with the 
winding-up of vertical field, by the weak differential rota¬ 

tion dz n of our disk model. While the inner domain boundary 
also shows an odd field symmetry, intermediate sections of 
the disk show a mixed parity, reflecting the presence of MRI 
modes with both even and odd vertical wave numbers. Quite 
remarkably, the MRI channels extend over several AU in the 
radial direction, although we remark that mis may be overly 
pronounced in the axisymmetric case, where the gr owth of 
parasitic mode s is likely to be artificially restricted ( jPessahj 
|& Chan||2008| ). Because we use a net-vertical field, the lin- 
ear MRI modes appear more pron ounced than in the other- 
wise comparable 3D simulations ofjDzyurkevich et al.jpOIOjl 
without a net Bz field. 

The poloidal velocity field plotted as black arrows is mostly 
disordered in the mrbulent regions but also shows some level 
of coherence in the upper layers, where a wind topology can 
be seen. While the wind is intermittent rather than steady and 
laminar, a general tendency for outflow is seen. This is con¬ 
sistent with similar observations of disk winds b eing launched 
from fully-MRI-active accretion d isks (see, e.g., Suzuki & In-j 
|utsuka|2009HPromang et al.|2013| l. We quantify me mass loss 
via the vertical disk surfaces by evaluating 


Afwind = StT 


pvg r sin 9 dr 


62 


t=ei 


(7) 


at the upper and lower disk surface, 9 = 0i, 02, respectively. 
Por model 0-b6, we find a net mass loss rate of Mwind = 
1.47 ± 0.37 X 10“® Mq yr-* (also cf. Table gbelow). 

We observe the formation of field arcs andmaterial flow¬ 
ing downward along field lines towards the arc foot points. 
This can be seen in the lower disk half around r — 2.8 au, 
and r = 4.3 au, for example, and illustrates the potential 
role of buoyancy instability in regulating the disk’s magneti¬ 
zation. Amplification of the magnetic field through the growth 
of channel modes, and their break up through the action of 
parasites, apparently leads to the formation of these locally 
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Figure 2. Visualization of the global disk structure for model 0-b6 with only Ohmic resistivity and /3p0 = 10®. The azimuthal fields reach peak strengths of 
~ 2.5 G in the active disk layers. Vertical outflow is intermittent, but regions of an ordered wind geometry can be identified (see upper half around 3.5 — 4 au). 



uprising field structures. As an alternative explanation, we re¬ 
mark that the dynamic evolution of such magnetic loops has 
previously been attributed to the effect of differential rotation 
onto the magnetic footpoints ([Romanova et al.|1998)l. 

In contrast, within the magnetically decoupled midplane 
layer, the magnetic field remains predominantly vertical, re¬ 
taining the initial field configuration. Near the dead-zone 
edge, the effect of Ohmic diffusion gradually declines. There, 
the azimuthal magnetic field is allowed to change its value, 
resulting in localized current layers. This is very similar to 
the situation obtained in the disks that include ambipolar dif¬ 
fusion, as discussed later. 

We illustrate the vertical disk structure in Fig. where 
we plot, in the upper panel, the total Reynolds and Maxwell 
stresses relative to the initial midplane gas pressure, pq. The 
profiles show the typical signa ture of a laminar dead zone and 
MRI-turbulent surface layers ([Fleming & Stone||2003[ [Oishi 
[et al. [[20071 [Gressel et al.[[2011[ [2012[ l, where the mechani- 
cal and magnetic shear-stresses lead to transport of angular 
momentum. Because of the relatively weak net-vertical mag¬ 
netic flux, the level of turbulence is modest, and the vertically- 
integrated dimensionless Maxwell stress is ~ 8 x 10“^. With 
the “viscous” mass accretion rate (i.e., that effected by inter¬ 
nal stresses) approximated by 



we estimate a value of Mvisc — 0.3 x 10 ® Mq yr at 
Tref = 2.5 au, if mass accretion were effected by turbulent 



- 8 - 6 - 4-20 2 4 6 8 



height z [H] 

Figure 3. Vertical profiles averaged over r S [2, 3] au and t = 25 orbits 
for model 0-b6 with Ohmic resistivity only and a midplane /3p0 = 10®. 
Negative values of the stresses (upper panel) are represented by thin lines. 


transport of angular momentum alone. This can be compared 
to the actual mass accretion rate, which we can easily com- 
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pute in our global model via 


-^accr = 27r / pv^ r sin 0 d0 
Je=0i 


(9) 


Applying a radial average over r G [2,3] au, and estimat¬ 
ing fluctuations occurring within a time interval spanning t G 
[75,100] yr, we find M^ccr = (0.14± 1.53) x 10"® M© yr"! 
(also see the last column of Table 1^. Clearly, this diagnos¬ 
tic turns out to be subject to strong fluctuations for model O- 
b6 - presumably due to the presence of strong radial motion 
within the MRI channel modes. Equation (j^ will however 
prove useful when evaluating the radial drift of material in the 
localized current layers seen in the AD-dominated disk mod¬ 
els. We note at this point that a simple energy conservation 
argument prevents the steady-state mass loss rate (to infin¬ 
ity) through a magneto-centrifugal wind being larger than the 
accretion rate through the disk. The large value of M^ind re¬ 
ported above relative to Maccr suggests that the win d loss rate 
is not yet converged. Indeed we note that |BS13| report that 
increasing the vertical sizes of their shearing boxes leads to a 
reduction of the wind mass loss ra tes. A similar conclusion 
is reached by|Fromang et al.j(|2013|) for winds launched from 
turbulent disks. It appears that obtaining accurate and phys¬ 
ically meaningful estimates for the mass fluxes through the 
upper boundaries of the disk will require simulations that are 
truly global in the vertical direction, such that the fast magne- 
tosonic point is contained within the simulation domain rather 
than outside of the boundary as it is at present for all of the 
models that we present in this paper (ensuring that the wind 
launching region is causally disconnected from the imposed 
boundary conditions). 

In the lower panel of Fig. the gas pressure, kinetic en¬ 
ergy and magnetic pressure, are shown. Within the dead-zone 
layer, between z = ±4 iJ, the disk is supported by gas pres¬ 
sure alone, which follows a Gaussian profile. Owing to the ad¬ 
ditional magnetic pressure support within the active layer, the 
effective scale height increases roughly twofold t here. N ote 
that unlike seen in the isothermal simulations of |BS13| cf. 
their fig. 3, the magnetic pressure does not significantly ex¬ 
ceed the value of the gas pressure in large parts of the disk 
corona. We attribute thi s difference to the dissipation heating 
( [Hirose & Turner|2011| l present in our simulations, which we 
note, however, may be overly pronounced in the axisymmet- 
ric models. The kinetic energy only becomes important very 
close to the vertical boundaries of our model. This is also re¬ 
flected in the position of the Alfven point, which is relatively 
poorly constrained and which we infer as (7.60 ± 0.45) iJ. 
Values for these vertical points are listed in Table for all 
models. We report time- and space-averaged values obtained 
for r G [1, 5] au, and we have appropriately taken into account 
the flaring of the disk. Because of the turbulent character of 
the upper disk layers, we cannot obtain a good estimate for 
the location of the base of the wind in model 0-b6. 

3.3. Laminar disk models with surface winds 

We now discuss models that include both O hmic re sistiv- 
ity and ambipolar diffusion. A key finding of BS13| is that 
the inclusion of ambipolar diffusion inhibits linear growth of 
the MRI in the intermediate disk layers, that is, the regions 
that are not affected by Ohmic resistivity. Since the ambipo¬ 
lar diffusion coefficient depends on the plasma parameter, one 
might expect that the effect of AD is less severe for high val¬ 




height z [H] 

Figure 4. Same as Fig.I^ but for model OA-b5 including both Ohmic resis¬ 
tivity and ambipolar diffusion, and with a stronger net-vertical field. Note the 
different scale of the abscissa in the upper panel, reflecting the reduced level 
of “viscous” transport. Dotted lines indicate the wind base at z = di^b- 


ues of /3p. To test this, we have run simulations OA-b6, and 
OA-b7, with /3po = 10®, and 10^, respectively. 

3.3.1. Ambipolar dijfusion with weak vertical fields 


In accordance with the /3po = cx) simulations of BSl3 we 
And very low levels of turbulence in these simulations (ct. Ta¬ 
ble 1^. As a consequence, the radial mass transport via ac¬ 
cretion (see Maccr values) is found to be negligible. At the 
same time, because of the weak vertical held, the magneto¬ 
centrifugal wind is absent in model OA-b7, and only poorly 
developed (and intermittent) in model OA-b6. We conclude 
that the corresponding held amplitude of about 10 mG (at 
1 au) can be assumed as a minimal requirement for signif¬ 
icant mass transport by magnetic effects. In the following, 
we accordingly focus on models derived from our fiducial run 
OA-b5, with /3po = 10®, that is, B,o = 31.5mG (4.2mG) at 
1 au (5 au). 

3.3.2. The fiducial model 

In Figure|^ we show vertical profiles of the Rf stress com¬ 
ponents for our fiducial model, OA-b5, where owing to the 
stronger net vertical held the MRI is ultimately suppressed 
by ambipolar diffusion and a laminar wind solution is estab¬ 
lished instead. Dashed vertical line s indicate the base of the 
wind, which we define according to|Wardle & K6nigl|(|1993|l 
as the vertical position, z = Zb, where the rotation becomes 
super-Keplerian. Since we define our azimuthal velocity, 
as the deviation from the Keplerian background flow, this 
can be verified in Fig. [^by the fact that the Reynolds stress, 
= pvp v^, vanishes at this point. The same is true for 


the vertical tensor component, 
nient to define the wind stress, 


which makes it conve- 


T i _ ry-iMaXW 


^ Z(f} 


_ /yiMaxw 
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at this point. For the i?(/) component of the stress tensor, re¬ 
sponsible for redistributing angular momentum radially, we 
furthermore obtain average values within the disk body z G 

[-Zb, +Zh], 


Tr4, = 


1 

‘2zh J- 


r 


Maxw 

Rep 


+ r^7")dz, (11) 


which we list separately in Table where values have been 
normalized to the midplane pressure, pq. Already for a mid¬ 
plane plasma parameter as low as 10^, the w ind stre ss exceeds 
the Ref) (accretion) stress by a factor of four. BS13 also report 
that the wind stress exceeds the accretion stress; for their fidu¬ 
cial run, the discrepancy is bigger than an order of magnitude. 

In the lower panel of Figure|^ we see that in the absence of 
MRI turbulence, the hydrostatic balance extends further up in 
the disk and the gas pressure remains the dominant stabiliz¬ 
ing force all the way up to the base of the wind. Even in the 
magnetically dominated wind layer, magnetic pressure gradi¬ 
ents remain moderate and the scale height of the disk remains 
roughly constant up to Zb. 

In Fig.l^ we show a zoom-in of the inner part of the global 
field topology in our fiducial simulation OA-b5. The up¬ 
per panel shows the resulting field configuration early on, 
when resistively modified MRI eigenmodes appearing inside 
r = 1.5 au are still clearly discernible. Localized channel 
solutions form at z « 3 H near the interface between the 
magnetically decoupled region, characterized by Aq ^ 1 
and purely vertical field lines with « 0, and the AD- 
dominated intermediate layer. In this small region, both El- 
sasser numbers are above unity (cf. Eig.[^, allowing for linear 
growth of MRI channel modes. This transition layer is charac¬ 
terized by a zigzag-shaped abrupt change in the field lines (see 
Eig. I^below), associated with strong current sheets. At this 
point, the laminar wind solution is already well established in 
the EUV ionized surface layer, and has spread throughout the 
radial extent of the simulation domain. 

While the disk wind itself has already reached a steady 
state at this point in time, the strong field layers continue to 
evolve. In the lower panel of Eig. H we accordingly show a 
later evolution stage, where the relicdVlRI channels have mor¬ 
phed into strong azimuthal field belts. These field belts can 
be understood as “undead zones”, that is, a region that does 
not sustain MRI, but where diffusion can be balanced by the 
winding-up of preexis ting radial field via differential rotation 
( [Turner & Sano|2008| . It appears that the initial amplification 
of the magnetic held through the growth of channel modes 
causes the ambipolar diffusion coefficient to increase until it 
quenches further development of the MRI. This notion is sup¬ 
ported by a reference simulation, where the diffusion coeffi¬ 
cients Tjo and ryAD were \\t\d fixed at their initial value, and 
where the MRI continues to produce quasi-turbulent fields 
akin to the ones seen in Eig. for the Ohmic-only case. The 
reversal of the field direction seen at r ~ 1.5 au is a relic of 
the local nonlinear development of the MRI modes early on 
(as may be seen in the upper panel at r ~ 0.75 au). This 
interface between the azimuthal field belts of opposite parity 
moves radially outward at a speed of « 0.01 auyr“^. 

The emerging current layers are directly related to 
resistively-modified vertically-global MRI channel modes - 


Latter et al. ( 20I0| l, and Fig.|^in Section 3.5 With MRI 


channels potentially spanning large radial extents, as in the 
Ohmic run shown in Fig.|^ such layers may be fed diffusively 
with magnetic field from the MRI-active inner disk. This is 



0.5 1.0 1.5 2.0 



0.5 1.0 1.5 2.0 

radius R [au] 


Figure 5. Field topology of our fiducial simulation at different evolution 
times. The azimuthal magnetic field (color) has been restricted to values 
|Bi^| < 125 mG for clarity; peak values are a few hundred mG. We also 
show projected magnetic field lines (white) and velocity vectors (black). Ad¬ 
ditional lines indicate the position, of the wind base (dot-dash), and the 
radial location of the profiles plotted in Figs.[^and^(dashed lines). 


supported by our present models, where current layers appear 
to emerge from inner disk regions. What determines the ex¬ 
act shape of the surviving field configuration at late times is 
presently unclear. We speculate that the particular realization 
seen in the lower panel of Fig.|^is not necessarily a unique so¬ 
lution given the specific disk geometry and ionization model, 
but may to some extend depend on the early evolution history. 
However, simulations that were initiated with different ran¬ 
dom seeds, but were otherwise identical, only showed minor 
variations in the final field appearance. The configurations 
observed at late times are moreover found to remain quasi¬ 
stationary, at least during the evolution time of 100 yr that we 
have investigated. 
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Table 2 

Summary of simulation results. 



2b 

ZA 

mHcyn 

^R(f> 

'pMaxw 

-^R4> 

ryMaxw 

Z(f> 

34wind 

^accr 


[H] 

m 

[10-%] 

[10-%] 

[10-%] 

[lO-«M0yr-i] 

[lO-»A40yr-l] 

0-b6 

— 

7.60 ± 0.45 

6.87 ± 14.4 

7.44 ± 0.95 

— 

1.58 ±0.59 

0.04 ± 1.91 

OA-b5 

5.23 ±0.07 

7.31 ± 0.17 

3.63 ±0.19 

2.22 ±0.06 

9.82 ±0.08 

0.79 ±0.01 

0.43 ±0.01 

OA-b6 

7.22 ± 0.48 

6.94 ± 0.37 

-0.21 ± 0.18 

0.79 ±0.06 

0.58 ±0.06 

0.19 ±0.06 

0.09 ±0.02 

OA-b7 

7.31 ± 0.70 

6.39 ±0.16 

0.07 ± 0.13 

< 0.01 

0.22 ±0.01 

0.03 ±0.01 

0.00 ±0.03 

OA-b5-d4 

5.27 ± 0.07 

7.33 ±0.18 

0.11 ± 0.30 

2.88 ±0.11 

9.88 ±0.12 

0.76 ± 0.02 

0.34 ± 0.04 

OA-b5-flr 

4.81 ± 0.03 

6.90 ±0.31 

0.26 ±0.21 

1.78 ±0.02 

14.3 ±0.02 

1.57 ± 0.01 

0.64 ± 0.02 

OA-b5-flr-nx 

4.78 ± 0.03 

7.50 ±0.30 

2.28 ±9.24 

1.87 ± 0.04 

13.0 ±0.04 

1.58 ±0.01 

0.64 ± 0.03 

OA-b5-nx 

5.10 ±0.04 

7.34 ± 0.13 

0.94 ± 9.29 

1.89 ±0.06 

7.84 ± 0.02 

0.67 ± 0.01 

0.30 ±0.04 


The vertical position of the base of the wind, z^,, and the Alfven point, Zj\, are found independent on the radial location when measured in local scale 
heights, H. The viscous accretion stresses are vertically integrated within \z\ < z^, - note the different units for Reynolds and Maxwell stresses. 
The wind stress, is measured at 2 = ± 21 ,. All stresses depend weakly on radius; listed values are at r = 3 au. 


3.4. Typical wind solutions 

Perhaps the most noteworthy result from run OA-b5 is the 
spontaneous emergence of the expected “hourglass” geometry 
for the magnetic field, allowing magnetic tension forces to 
accelerate the wind. For the adopted vertical field polarity, the 
field must bend radially outward near the upper and lower disk 
surfaces, and the azimuthal field should point in the positive 
direction in the lower hemisphere and reverse direction above 
the midplane, as observed in the upper panel of Figure]^ 

3.4.1. Robustness of the emerging wind geometry 


The shearing box simulations of BS13 have reflectional 


symmetry in the radial direction, that is, with respect to mir¬ 
roring about a; = 0, and hence contain no information about 
the location of the star. Instead of giving rise to a physical 
wind solution, their fiducial run produces a radial field that 
possesses odd-symmetry about the midplane, with an inward 
pointing field in one hemisphere and an outward pointing one 
in the other. To test the robustness of the emerging solution 
in our setup, we initiate several instances of run OA-b5, us¬ 
ing different random number seeds when perturbing the initial 
velocity field, and in each case we recover the correct wind so¬ 
lution with outward bending field lines. Our initial model has 
a net-vertical field with a radial dependence such that the mid¬ 
plane value of /3p remains constant. Because the gas pressure 
itself decreases with radius, the corresponding Bz{R) (with 
QrBz < 0) results in an outward magnetic pressure force, 
which moreover has a stronger acceleration effect in the low- 
density upper disk layers. While a vertical flux distribution 
that has the flux centrally concentrated can be expected when 
accounting for inward advection and outward diffusion of flux 
( Okuzumi et al.|[2014| l, our particular choice is of course ar- 
bitrary. Radial magnetic gradients are furthermore known to 
affect the outflow’s mass loading in axis ymmetric calculatio ns 
where the disk is a boundary condition ( Pudritz et al.|2006) l. 

In our simulations, the pressure gradient related to B^ (R) 
may potentially affect the direction toward which the field 
lines first bend from their starting configuration. We have 
tested this by running reference models with different mag¬ 
netic pressure gradients. If dnB^ > 0, the initial condition 
has an unbalanced inward pressure force, and we do indeed 
find the field lines to spontaneously bend inward', this hap¬ 
pens simultaneously in the upper and lower hemisphere of the 
disk, such that the overall symmetry that is obtained is again 
even. The launching of the inward wind is restricted to the 


inner disk, probably because the vertical field lines become 
too stiff (owing to increasing with radius) to be suitable 
for wind launching in the outer disk. 

For a neutral situation with dnB^ = 0, we still observe out¬ 
ward bending of the field lines. Starting from the inner radial 
domain, the outward propagating establishment of the wind 
region is found to stall at some radius, whereas the wind was 
quickly established throughout the entire domain in model 
OA-b5. As before (for the case of an outward magnetic pres¬ 
sure gradient), this is probably related to the field lines be¬ 
coming too strong to support the wind launching mechanism 
at large radii. The wind indeed stalls further out in a run 
with weaker overall net-vertical field. This poses the ques¬ 
tion whether the vertical profiles of Ao(^), and Aad( 2) that 
we obtain from our ionization model put strong constraints 
on permissible vertical field amplitudes as a requirement for 
wind launching. 

A possible reason for the outward bending, even in the case 
of neutral magnetic pressure forces, may be that the vertical 
shear present in the global models (because of the baroclinic 
conditions) naturally bends the azimuthal component of the 
field in the correct direction required by the physical wind 
solution. We have checked that the outward bending of the 
vertical field lines is however equally seen in (strongly flar¬ 
ing) globally isothermal models that do not have any vertical 
gradients in the rotation velocity or any radial gradients in the 
vertical magnetic field. The outward pointing configuration 
is energetically favored because setting it up involves diluting 
rather than concentrating magnetic flux. 

Our simulations demonstrate in any case that global models 
spontaneously develop the correct field geometry, although 
we caution that this conclusion needs to be tested in future 
simulations that moreover adopt improved boundary condi¬ 
tions for the magnetic field at the upper and lower (and po¬ 
tentially the inner radial) surfaces of the simulation domain. 
Ultimately, the emergence of the wind geometry will have to 
be studied in simulations that do not start from well-motivated 
(but nevertheless arbitrary) initial conditions, but do account 
for th e formation of the PPD from a collapsing molecular 
cloud ( |Li et al.|2014[ ). 

Lastly, for the parameters that we have considered here, we 
do not fi nd a self-limiting of th e wind via shielding of FUV 
photons ( Bans & Konigl 2012 1 . Th i s is co nsistent with the 
T Tauri scenario of|Panoglou et al.| (|2012|l, who performed 
disk wind chemical modeling for vmous protostellar evolu- 
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wind solution at R=1.0au 




height z [H] 

Figure 6. Vertical profiles of the velocity components (upper panel) and the 
magnetic field components (lower panel) projected to cylindrical components 
for run OA-b5 at spherical radius r = 1 au (see the inner dashed line in 
Fig. pp). The Alfvenic point with respect to uaz = is marked by 

filleocircles. Dotted vertical lines indicate t he base of the wind, defi ned as 
the point where becomes super-Keplerian jWardle & K6nigl|1993^. 


tion stages. The authors hnd the shielding of FUV photons 
to be important for the Class I and Class 0 cases. Their 
T Tauri case, with mass flow rates comparable to our hdu- 
cial model, is sufficiently FUV-ionized to reduce ambipolar 
diffusion enough so that the neutrals are swept up in the wind 
out to at least a radius of 9 au. 


3.4.2. Comparison with local models 

Except for the reversal of the horizontal magnetic field com¬ 
ponents (which at this stage appears to occur rather gently 
instead of within a narrow current layer), the field structure 
observed i n the u pper panel of Fig. is similar to that de¬ 
scribed by|BS13|and shown in their figure 10. Near the mid¬ 
plane, where the dominance of the Ohmic resistivity retards 
the growth of currents, the field is dominated by the vertical 
component. Moving to higher altitudes, where the magnetic 
coupling increases, the field bends outwards because a radial 
velocity is generated through the azimuthal force balance be¬ 
tween the Coriolis force and magnetic tension. We enter the 
FUV layer at disk altitudes \z\ > A.5H, coinciding with the 
base of the wind and the wind itself, as illustrated by the ve¬ 
locity vectors. 

At later times, the lower panel of Fig.l^shows that the vary¬ 
ing polarity of the strong azimuthal helclbelts gives rise to two 
distinct field line geometries. In the inner region (r ^ 1.5 au) 
the orientation of the held belts is the same as the one in the 
wind, and the vertical held lines bend smoothly at the transi¬ 


wind solution at R = 2.0 au 




height z [H] 

Figure 7. Same as Fig.|^but for r = 2 au, where the horizontal field belts 
have opposite parity compared to the field within the wind layer. Inlays mag¬ 
nify the velocity profiles near the strong current layers, which reach a fair 
fraction of the local isothermal sound speed, Cg = 1.1 kms“^. 


tion into the wind base. Because some negative has dif¬ 
fused into the Ohmic-dominated midplane, there is a weak 
current layer forming at z ~ —2.15 H. In contrast, further out 
at r 1.5 au, the horizontal-held belts are anti-aligned with 
the magnetic held direction of the wind. This can also be seen 
in the curvature of the held lines, which are concave (towards 
the star) in this region. To connect the different layers, the az¬ 
imuthal and radial held components have to change their di¬ 
rection within two thin current layers located at z ~ ±3.2 H. 
To study these conhgurations in more detail, we plot sepa¬ 
rate vertical prohles for r = 1 au, and r = 2 au in Figs. 
and It] respectively. For easy reference, the radial positions 
are fmthermore indicated in Fig. by means of dashed lines. 
The isothermal sound speeds, for the initial disk model, are 
Cg = 1.5 km s“^, and 1.1 km s“^ at the two respective radii. 

The wind prohle in the inner disk (see Fig. is qualita¬ 
tively similar to the wind solutions previously found in the 
quasi -ID si mulations with “even-z” symmetry, plotted in hg. 
10 of BS13| The current layer is much weaker in our case, and 
the dip in is barely visible in the left inset of Fig. This 
may be related to differences in the ionization model^esult- 
ing in Aad > 1 in the region around |z| ~ 3iT, in our disk 
model. Differences in the particular diffusivity prohles also 
explain the pronounced peaks in in this “unde ad” tra nsi- 
tion layer, which are not seen in the simulations of |BS13| 

The prohle at r = 2 au, with double held reversal in B^, 
is illustrated in Fig. Because of the strong current layers 
that are emerging in this case, this conhguration resembles 
more closely the situation seen in the local models (again cf. 
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Figure 8. Vertical structure of diffusively modified MRI channel modes, with 
Ohmic resistivity only. The solution that has been extracted from the non¬ 
linear global simulation (light color) agrees markedly well with the linear 
eigenmode solution (dark color). 


fig. 10 in|BS13|l. Unlike in their simulation, we however ob¬ 
serve two current layers, one on each side of the disk. These 
arise because the development of MRI channel modes, and 
their subsequent evolution into “undead” azimuthal field belts 
with opposite polarity to the background wind, forces the hor¬ 
izontal field to change direction three times. Strong current 
layers occur at altitudes where the field and gas couple more 
strongly, namely at \z\ ^ 3H. Although the occurrence of 
strong azimuthal field belts is not precluded in local simula¬ 
tions (as their existence may depend on the details of the ion¬ 
ization model), it seems likely that the alternating field belts 
we observe will only be a robust feature of global simula¬ 
tions. Apart from their multiplicity, the cur rent layers appear 
very similar to the ones presented in|BS13| The radial veloc¬ 
ity Vr shows a dip, while is modulated as a sine (the slight 
offset of in our simulation is due to the radial pressure 
support). The characteristic modulation (reflecting the abrupt 
kink in the held lines) is clearly seen in and B^, which 


agree very closely with the B^, and By obtained by BS13 
The emergence of such strong current layers naturally begs 
the question of their stability. Before we however discuss this 
issue in section [T6l we first want to take a look at the current 
sheets’ likely origin as remnants of MRI eigenmodes. 


tization is lowest and the MRI wavenumbers are the highest, 
explaining the similarity of the solutions for constant rj, and 
spatially varying rj = ri{z). 

Here, instead of a s tratified shearing bo x framework as pre¬ 
viously employed b y [Latter et al. (|2010| l , we u sed the frame¬ 
work introduced by McNally & Pessah ( |2014| ) to capture the 
full vertical structure of the baroclinic initial condition. The 


eigenmode analysis uses the linear equations described in Mc¬ 


Nally & Pessah (2014 Section 5) with the addition of terms 
for spatially varying Ohmic resistivity, and using the func¬ 
tional forms and coefficients appropriate for the c ylindrical 
temperature stru cture of the disk initial condition ( McNallyl 
& Pessah|2014[ Section 3.3). The equations were discretized 


in real space with Chebyshev cardinal functions on a Gauss- 
Lobatto grid with 1000 points with the same vertical extent 
as the global simulation (±8 scale heights). Boundary con- 
ditions were enforced by the “boundary bordering” method 
( Boyd|2000 1 forcing the magnetic held perturbations to zero 
at the vertical edges of the domain. The analysis was per¬ 
formed at radial position i? = 1 au where the Ohmic resistiv¬ 
ity profile, rjo{z), used was extracted from the simulation ini¬ 
tial condition and interpolated to the grid used for the eigen¬ 
mode calculation. 

We moreover obtained eigenmode solutions that addition¬ 
ally included a fixed pad ( 2 :) profile, and these showed a 
nearly indistinguishable mode structure. Because of the ex¬ 
tra non-linearity introduced by the /3p dependence of pad( 2), 
we however restrict our comparison to the Ohmic-only case. 
In Figure we compare the early stage of our fully non¬ 
linear global simulation with the fastest growing eigenmode. 
The localized held reversals are also clearly seen in the upper 
panel of Fig. which shows the early evolution of model 
OA-b5 (however including AD). To isolate the linear MRI 
mode from the overlaid wind solution, we subtract respec¬ 
tive held profiles from two distinct snapshots during the linear 
growth phase of the MRI. This is possible because the wind 
configuration emerges quickly, and then provides an approxi¬ 
mately stationary background in which the MRI grows. Dif¬ 
ferencing two snapshots in time hence allows us to remove the 
wind part and extract the exponentially-growing MRI eigen¬ 
mode. While growing perturbations near the midplane are 
high-wavenumber and thus efficiently damped, the relative 
held strength in the disk halo is too strong to allow for the 
MRI modes to fit within the available domain. This only 
leaves an intermediate disk region near j^l ~ 3.5 iT for MRI 
channel modes. In the Ohmic-tAD simulations, this corre¬ 
sponds to the position of the “undead” held belts (see lower 
panel of Fig. with adjacent current layers. In those simu¬ 
lations, the back-reaction of the growing eigenmodes on the 
Pad ( 2 ) profile makes the MRI self-quenching, leaving the 
current layers as a remnant of the initial MRI channel mode. 


3.5. Relation between current layers and channel modes 

To draw a closer connection between localized MRI chan¬ 
nel modes and the current layers observed in our simulations, 
we compare the structure of the fastest growing MRI eigen¬ 
mode with the initial evolution of our standard setup with 
/3p0 = 10®. The mode structure we And is qualitatively sim¬ 
ilar to that shown for a spatially constant O hmic resistivity 
(cf. the appendix of jFromang et al.||2013]l, and rese mbles 
the eigenmodes shown in figure 4 ofjLesur et al.|(|2014|l, who 
performed Hall-MHD calculations that also included variable 
Ohmic resistivity. Dissipation modifies the MRI eigenmodes 
predominantly near the midplane, where the relative magne¬ 


3.6. Stability of current layers and horizontal field belts 

While we And the current layers to remain stable and long- 
lived in model OA-b5, some of the other simulations which 
we will discuss below also show signs of instability, leading 
to less regular fields. 

3.6.1. Effect of disk flaring 

An example of this is illustrated in Fig. where we show 
the vertical velocity, (upper panel), and tne azimuthal mag¬ 
netic held, B^ (lower panel), at the end of the flaring disk 
simulation OA-b5-flr. Both color scales have been restricted 
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Figure 9. Close-up showing the thin current layer of the flaring disk model 
OA-b5-flr breaking up into a KH-type instability. Top: vertical velocity, . 
Bottom: azimuthal magnetic field strength, and poloidal magnetic field 
lines. The wavenumber of the perturbations is found to depend non-trivially 
on radius. Only for r ^ 1.5 an, the nonlinear stage is clearly identified. The 
KH pattern is furthermore imprinted onto the mass-loading of the wind. 


to highlight features in the current layer. In this model, the 
surface density interior to 1 au is smaller than in OA-b5, and 
beyond 1 au it is larger, altering the ionization structure suffi¬ 
ciently to put the current layer in a different regime. In partic¬ 
ular, the increased ionization in the inner disk apparently al¬ 
lows the azimuthal field belts that form at intermediate heights 
to reach larger amplitudes. This, combined with the fact that 
these field belts are of opposite polarity relative to the back¬ 
ground wind field, causes the two current layers located at 
z ~ ±3H to be stronger and to support narrow, radially di¬ 
rected flows, or “jets”, similar to those in run OA-b5, as illus¬ 
trated in Figure!^ 

Inside r ~ l.Tau, the velocity near the current layer shows 
a distinct radial wavenumber, suggesting that the shear associ¬ 
ated with the radial jets leads to growth of a Kelvin-Helmholtz 
(KH) instability. The pattern, however, remains more or less 
stationary, indicating that the instability is not growing signifi¬ 
cantly into the nonlinear stage. Although not obvious in Fig|^ 
wave-like perturbations to the magnetic held are also present 
in the current layer. The possibility of current-driven tearing 
modes being present also arises, but analysis of the shear rate 
associated with the radial jets suggests that the observed in¬ 
stability is driven by KH modes. The velocity associated with 
the shear is super-Alfvenic, and is therefore expected to stabi¬ 
lize the teari ng modes that may in principle grow in a shearing 
background ( |Chen et al.|l997| . 

The wavelength of the velocity perturbations and magnetic 
held changes abruptly for r 1.4 au, where now a character¬ 
istic vortex-pair pattern emerges. The exact reason for this 
abrupt transition is unclear, since most characteristic quanti¬ 
ties, such as the Elsasser numbers Aq, and Aad, only depend 
weakly and continuously on radius. The most likely explana¬ 


tion may be the radial variation of the azimuthal held (inside 
the “undead” belts), which drops in strength around this radial 
location. This is accompanied by a weakening and widening 
of the current layer, and a significant reduction in the inflow 
velocity of the radial accretion flow associated with the cur¬ 
rent layer, and a corresponding reduction in the shear relative 
to the background fluid. This suggests that the abrupt change 
in the nature of the disturbance observed at radius ~ 1.8 au is 
unlikely to be driven by a KH instability, but the weakening of 
the shear may allow a tearing-mode instability to develop in¬ 
stead. The appearance of disturbances that have the character 
of growing magnetic islands provides some support for this 
interpretation, but this is clearly a tentative explanation that 
will require further investigation in future work. The pres¬ 
ence of a background wind, Keplerian shear, ambipolar dif¬ 
fusion and a radially-directed inflow associated with the cur¬ 
rent layer make an analysis of this problem somewhat com¬ 
plicated, and beyond the scope of this paper. We note, how¬ 
ever, that the local stud ies on the development of par asitic 
modes in resistive MHD (|Latter et al.|2009[|Pessah|2010|l may 
provide some insight, eve n though ambi^ lar diffusion is ne¬ 
glected in those studies. Pessah ( 2010| l has analyzed MRI 
parasitic modes in the presence of finite Ohmic resistivity and 
viscosity. In this paper, Aq, is identified as the relevant di¬ 
mensionless number (see their figure 11), with tearing modes 
being the dominant parasitic modes for Aq ^ 1, and Kelvin- 
Helmholtz modes dominating otherwise. Considering both 
Aq, and Aad, we assess that our simulations are situated near 
this transition point. 

3.6.2. Ejfect of reduced dust fraction 


Figure 10 shows the final appearance of the entire disk for 
model OA^5-d4, which has a reduced dust fraction of 10“"^ 
compared to a value of 10“^ in model OA-b5. On one hand, 
as the reader may check from Table the two models are 
rather similar in their bulk properties, especially those related 
to the wind. This is unsurprising, since the reduced cross- 
section for recombination on grains affects the ionization frac¬ 
tion significantly only at large gas columns (see Fig.[2l, below 
the base of the wind. The grains are less important for the 
ionization state in the FUV layer, where the electrons are so 
abundant that many remain free once the grains charge to the 
Coulomb limit ( |0kuzumi|20091 l. 

In contrast to our flducial model, the transition between the 
base of the wind and the Ohmic-dominated disk interior is 
very chaotic in model OA-b5-d4, whereas it is much more 
stable in model OA-b5. It appears that the larger values of 
Aq and Aad in this model at heights z ^ ±2.5iT allow the 
channel modes to grow further into the nonlinear regime, such 
that parasitic modes are able to develop and cause the held 
belts to break up into regions with locally coherent held. The 
held amplification during this phase means ambipolar diffu¬ 
sion strengthens, as in run OA-b5. Eventually the local evolu¬ 
tion is driven by a combination of field winding and diffusion. 
The end result is the creation of local “undead” patches that 
evolve on timescales that are comparable to the run duration. 

The different evolution of the horizontal field belts observed 
in runs OA-b5-flr and OA-b5-d4 compared with OA-b5 show 
that the evolution of the disk interior depends critically on the 
precise disk model that determines the Elsasser number pro¬ 
files. In particular, the ionization state in the region between 
the Ohmic-dominated disk midplane and the AD-dominated 
intermediate disk layers depends sensitively on various model 
parameters. This suggests that a careful parameter study is re- 
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Figure 10. Same as figure Fig.^ but for model OA-b5-d4 with reduced dust fraction of 10“^, and showing the whole disk. Residual MRI fields form irregular 
patches below the wind layer leading to spatially intermittent current sheets. 


quired to obtain a reliable picture of how PPDs evolve over 
their life-times - especially since their properties change sub¬ 
stantially over these time scales. Such a parameter study 
would appear to be most warranted when additional physics 
is included in the model, such as the Hall effect and more 
realistic thermodynamics. In a situation where the gas is al¬ 
lowed to be heated by resistive effects, the ionization balance 
may be shifted towards higher conductivity in the vicinity of 
the current layers, potentially leading to intermittent behavior 




E aure, Fromang, & Latter|2014| McNally, Hubbard, Yang, & 

ac Low|2014r --- 


3.7. High-wavenumber instability at FUV transition 

A distinct instability emerges in the magnetically domi¬ 
nated corona of the laminar disk (see Fig.[TT]l. The instability 
is found to grow close to the smallest available wavenum¬ 
ber, and is only present in the highly-resolved axisymmet- 
ric simulations. The perturbations are roughly aligned with 
the outward-pointing field lines of the laminar wind solu¬ 
tion. Similar patterns have been observed in 2D-axisymmetric 
sh earing box si mulations with strong vertical fields performed 
by |Moll| ( |2012[ ), who also finds field-aligned structures on the 
shortest length scales resolved on the grid. The ins tability 
also bears resemblance with the “striped wind” seen by|Lesur, 
Ferreira, & Ogilvie p013| l. Similar to the instability seen by 
Moll| ( ^12| l, we find a rnoderate modulation of the vertical 
mass flux by the instability pattern. We however do not find 
material falling down along the field lines as reported there 
for certain parameters. 

We identify three candidate mechanisms for creating such 
disturbanc es: (i) the vertical-shear instabi l ity inv estigated re¬ 
cently by Nelson, Gressel, & Umurhan (|2013[), potentially 



radius R [ou] 

Figure 11. Color-coded azimuthal component of the flow vorticity (at late 
time in model OA-b5), tracing high-frequency perturbations in the magneti¬ 
cally dominated disk corona. Perturbations originate at the base of the FUV 
layer (dashed line), and are aligned with the field lines of the wind solution. 


modified by anisotropy from field-line tension, (ii) a Kelvin- 
Helmholtz type instability related to the vertical gradient of 
the horizontal velocity near the base of the wind, and (iii) a 
buoyant interchange instability of the azimuthal field. It is 
important to note that the KH-type instability can create the 
elongated aspect ratio seen in Fig. El when perturbations are 
catried-away by the overlaid laminar wind velocity - and we 
indeed find this to be the case in situations where the instabil¬ 
ity occurred before the wind was established. The conclusive 
hint that a Kelvin-Helmholtz instability causes the pattern in 
Fig. [TT] is a strong velocity gradient at the exact location of 
the transition into the FUV layer. This gradient can be clearly 
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identified in Fig. [^(situated just above the wind base at Zb). 
The shear layer manifests itself as an additional steepening of 
the wi nd prof iles, and can also be seen in the corresponding 
plot of BS13 (see their figure 10). 

The physical origin of this sharp velocity gradient can be 
traced back to a peak in the ambipolar part of the electromo¬ 
tive force in the induction equation. Because the coefficient 
?7AD(f, 0) sharply drops at the position of the FUV interface, 
the term V X (pad [(VxB) X b ] X b) can become significant 
in the induction equation even for smoothly varying (or con¬ 
stant) B. In the steady state, this curl of the ambipolar elec¬ 
tromotive force needs to be balanced by the induction term, 
V X (v X B), requiring as a consequence a significant shear 
gradient in the horizontal velocity (again assuming a smooth 
B). Given its strong dynamical effect, this raises the ques¬ 
tion whether the coincidence of the wind base and the FUV 


interface is in fact accidental as claimed by BS13 
In our disk model, the FUV layer is imp osed ad hoc, with 
a prescribed column density (Section The resulting in¬ 

terface is smeared out to a prescribed width. Reflecting the 
predominant direct origin of the FUV photons, while scatter¬ 
ing is expected to remain negligible, the tapering is done in the 
radial direction, and as a function of column density (rather 
than position). As a consequence, the resulting pad ij-, 6) may 
develop sharp steps on the numerical grid in the 0 direction - 
in turn resulting in sharp velocity gradients. We accordingly 
ran a reference model, where the FUV transition was softer. 
In this model, we did not observe a strong velocity gradient 
near the wind base, and also found the instability to be largely 
suppressed. Also the flaring-disk model, OA-b5-flr, did not 
show any signs of the instability. This can potentially be ex¬ 
plained by the FUV transition following a concave line in this 
case, producing a different aliasing of the tapering function 
on the spherical-polar mesh. We conclude that the instabil¬ 
ity seen in Fig.f^is likely caused by the particular numerical 
parametrizationOT the FUV layer in our mod el, and is ent irely 
unrelated to t he “st riping” instability seen in Moll ( 2012| l and 
|Lesur et al.| ( 2013^ The instability may be be triggered by 
numerical effects at the grid scale, in a manner analogous to 
how physical secondary Kelvin-Helmholtz instabilities can be 
triggered by p urely numerical effects as a shear layer thins to 
the grid scale ( McNally et al.||20T2 i. We remark that in the 
case of this wind instability too the mechanism behind the in¬ 
stability is genuinely physical, and may occur in real PPDs 
where there is a sharp local change in the ionization state. 


3.8. Non-axisymmetric simulations 

The presence of secondary instability mandates investigat¬ 
ing the non-axisymmetric evolution. This is done exemplar- 
ily for the two models OA-b5, and OA-b5-flr, whose 3D 
counterparts OA-b5-nx, and OA-b5-flr-nx are also listed in 
Table showing only minor deviations from their respec¬ 
tive 2D runs. Given the moderate magnetic Reynolds num¬ 
bers in these regions, it is not per se obvious whether non- 
axisymmetric instabilities can transfer significant amounts of 
energy into modes with higher azimuthal wavenumber. To ad¬ 
dress this question, in Fig. we show vertically-averaged 
azimuthal power spectra orme magnetic field. The color¬ 
coding indicates the logarithm of the power spectrum of Bg 
as function of wavenumber, m, and radial position in the disk. 
In the upper panel, we show the result for model OA-b5-flr. 
With the exception of the Kelvin-Helmholtz features (between 
R = 1.5 — 2.0 au) also seen in Fig.|^ almost all power resides 
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Figure 12. Vertically-averaged azimuthal power spectra of the Bg magnetic 
field. Top: Model OA-b5-flr with a laminar wind. Bottom: A run w ith MRI 
active surface layers in a domain of slightly different radial size, from|Gressel| 
|et al.H2013J. Owing to the restricted azimuthal extent in both runs, the only 
wavenumbers present are multiples of four. 


in TO = 0. Even in the regions of the secondary instabili¬ 
ties, the power in modes with to > 0 remains very small. 
The power seen adjacent to the inner radial domain boundary 
is caused by reducing the diffusion coefficients there, which 
is a n att empt to mimic the MRI-active inner disk (see Sec¬ 
tion for details). 

This can be contrasted to the case of the classic l ayered 
PPD with MR I-active surface layers (see model Ml in Gres-| 
sel et al. 2013[), which we plot in the lower panel of Figure) 127 


Even with the laminar midplane layer, a significant amountof 
power is found for high-TO modes (note the different color bar 
in this figure). The lack of power near the radial boundaries 
is owing to a dissipative buffer region. |Gressel et al.| P013 i 
focused on the effects of a planet, which they inserted after 
the time shown here. We hence conclude that in the simu¬ 
lations that are dominated by axisymmetric configurations of 
AD-laminar winds, secondary instabilities do not lead to sig¬ 
nificant power in non-axisymmetric perturbations. This may 
potentially change in the presence of stronger current layers 
that are prone to nonlinear stage of parasitic modes as, for in¬ 
stance, can be seen for model OA-b5-d4 with a stronger dust 
depletion. When combined with stronger net-vertical fields 
and flaring, these dust-depleted disks may develop intermedi¬ 
ate layers that are on the verge of becoming turbulent. 


4. SUMMARY 

We have presented the first quasi-global MHD simulations 
of PPDs that include Ohmic resistivity and ambipolar dif¬ 
fusion. In agreement with the shearing box simulations of 
BSI3| we find that the inclusion of ambipolar diffusion has 
a dramatic effect on the evolution. Our main results may be 
summarized as follows. 
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1. In simulations where we include a weak vertical magnetic 
field and Ohmic resistivity, but neglect ambipolar diffu¬ 
sion, we obtain evolution that is consistent with the tradi¬ 
tional dead zone picture: surface layers that sustain MRI- 
driven turbulence and a quasi-laminar midplane region. 

2. The inclusion of ambipolar diffusion in models with weak 
magnetic fields, where the midplane value of /?p > 10^, 
results in quenching of the surface layer turbulence. Such 
disks display very weak internal motions and negligible 
angular momentum transport. 

3. Our most important result arises when we increase the 
strength of the magnetic field so that = 10®. This 
leads to a disk which is entirely laminar between 1-5 au, 
and, despite the still relatively weak field, to the launch¬ 
ing of a magnetocentrifugal wind from the upper ion¬ 
ized layers, resulting in a significant wind stress being ex¬ 
erted on the disk that drives inward mass flow at a rate 
> 10-®Moyr-i. 

4. We have tested the robustness of the wind solution by 
running multiple simulations with modified initial condi¬ 
tions. Because of the outward pressure force from the 
initial Bz{R), and the vertical differential rotation, we 
generally obtain magnetic fields that point outwards in 
the radial direction in both hemispheres, as required for 
a physical wind solution that extracts angular momentum 
from the disk. This is with the exception of a model with 
SrBz > 0, where an inward pointing wind was found. 

5. The calculations treat the time-varying absorption of the 
FUV radiation from the star and its vicinity in the inter¬ 
vening material, yet there is no sign that absorption by the 
wind limits the ionization of the launching layer for mass 
loss rates in the T Tauri range. 

6. Strong horizontal field belts are formed as remnants of 
stagnated MRI modes at intermediate heights in our disk 
models, and are likely maintained through a balance be¬ 
tween field winding and ambipolar diffusion. These are 
found to have field polarities that may be aligned or anti¬ 
aligned with the background field associated with the 
wind, affecting the number of current layers that arise in 
the disk because of local field reversals. The belts polarity 
may depend sensitively on the early evolution of the MRI. 

7. The current layers drive substantial accretion flows that 
are narrowly localized in height. The associated shear can 
be strong enough to allow small-scale Kelvin-Helmholtz 
instabilities to grow. 

8. The global evolution of the horizontal field belts is found 
to depend on the Ohmic and ambipolar Elsasser numbers 
(and hence the ionization fraction). In our fiducial model, 
the belts alternate in polarity over large radial extents, but 
otherwise appear to be stable and only evolve on time 
scales that are long compared to the dynamical time scales. 

9. Reducing the dust grain abundance by an order of magni¬ 
tude, or allowing the disk to have a flaring structure, leads 
to more complex evolution in which the belts break-up into 
smaller-scale islands of locally coherent horizontal field 
polarity. We hypothesize that this evolution results from 
the increased ionization fraction that allows MRI channel 


modes to grow further into the nonlinear regime, where 
they are prone to break up via parasitic modes. 


Our global simu lations are most comparable to the shearing 
box simulations of |BS13[ and so it is worth briefly comparing 
our results with theirs. Our results for weakly magnetized 
disks (/3p > 10®) with ambipolar diffusion and Ohmic resis¬ 
tivity compare well with their results for disks with no net ver¬ 
tical magnetic field, that is, no sustained wind and very weak 
turbulence and angular momentum transport. This indicates 
that there is a clear lower limit for the net vertical magnetic 
field required to drive an accretion flow that is in agreement 
with observations (given our restricted assumptions about the 
disk microphys ics). O ur results are in complete qualitative 
agreement with|BS13|for higher magnetic field strength: disks 
are laminar and a magnetocentrifugal wind is launched from 
the highly ionized surface layers. The associated wind stress 
is sufficient to drive an accretion flow that matches typical 
rates measured to be reaching the surfaces of young stars. 

Unlike earlier studies of magnetocent rifugal wind launch¬ 
ing from deeper towards the midplane (Blandford & Paynej 


1982||Pudritz & Norman|1983||Wardle & K6nigl| 1993)1, there 

IS no requirement for a strong magnetic field to be present 
in the disk when the wind is launched from high altitudes, 
since the requirements for field rigidity are met for globally 
weak fields in these low density regions. The appearance of 
long-lived horizont al field belts is not observed in the shearing 
box simulations of |BS13[ but this is because of modest dif¬ 
ferences in the disk models and assumed ionization physics 
(their model has a somewhat larger surface density at 1 au, 
and we include a direct component to the disk illumination 
by stellar X-ray and FUV photons). The models agree qual¬ 
itatively on the appearance of strong current layers, but the 
stability of these could not be examined by |BS13| as their 
models in which the field geometry gave rise to these cur¬ 
rent layers used a quasi-ID approximation. This latter point 
illustrates the primary difference between the global and local 
shearing box simulations. The global simulations appear to 
produce and sustain the correct physical wind geometry spon¬ 
taneously, whereas the symmetries associated with the shear¬ 
ing box normally lead to a radial magnetic field that points 
inwards in one hemisphere and outwards in the other, such 
that the correct geometry needs to be enforced. The potential 
role of the rather arbitrary initial conditions, and the vertical 
boundaries in setting the magnetic field topology in global 
models should be investigated. Our choice to impose a verti¬ 
cal field is motivated by the typical “hourglass” m orphology 
(e.g. figures 2 and 3 in Kunz & Mouschovias|2010| and refer¬ 
ences ther ein) emerging also in z oom-in simiulations of PPD 
forma tion (|Nordlund et al.|2014| l, and observed in molecular 
cores (jGirart et al.|2006|l; even though the topologies of PPDs’ 
magnet fields have not yet been unambigu ously determined 
( [Stephens et al.|20T4 Davidson et al.|20l4 l. 


5. CONCLUSIONS 

The results presented in this paper, combined with those by 
BS13| indicate that the traditional picture of a PPD hosting a 
dead zone in the shielded disk interior, with turbulent surface 
layers driven by the MRI, no longer holds. When ambipolar 
diffusion is treated, the flow becomes laminar, and accretion is 
driven by the magnetic stresses produced in launching a wind 
from the disk surface. Adding the Hall effect appears to fur¬ 
ther modify the flow by allowing a significant Maxwell stress 
to develop near the disk midplane through the combined ac- 
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tion of field winding and Hall drift (|Kunz 2008 Kunz & Lesur 


|2013[ [Lesur et^|2014t |Bai|2014| l, but does not lead to tur 

bulence. We are then left with a picture of PPDs in which the 
innermost regions sustain vigorous MRI-induced turbulence 
because of collisional ionization of alkali metals, the interme¬ 
diate regions between ~ 0.5 — 20 au are laminar with angular 
momentum transport occurring through a disk wind and per¬ 
haps through Maxwell stresses near the midplane (an effect 
that only arises if n • B > 0 and when including the Hall 
term), and an outer region where the Hall effect and Oh mic 
resistivity are sub-domina nt and weak turbulence arises s 
jmon et aL|2013[[Bai|2015| l. 

The picture painted above has important implications for 
dust evolution and planet formation. In principle, dust coagu¬ 
lation and settling in the laminar region should occur rapidly, 
and in the absence of a source of stirring the loss of small 
grains should affect th e spectral energy distribution ( SEP). 
This was examined by Dullemond & Dominik (20051, who 
showed that this rapid evolution of the SED is not in agree¬ 
ment with observations, and for many years this important 
result has provided indirect evidence for the existence of tur¬ 
bulent stirring in the planet forming regions of protoplanetary 
disks. In the absence of a magnetic origin for local turbulence 
this suggests that either another source of turbulence exists. 


such as the vertical shear instability (e.g. Nelson, Gressel, & 
Umurhan]|2013 1, or that small grains are continuously deliv 


ered to the laminar region from the outer and/or inner turbu¬ 
lent regions. 

The results in this paper are consistent wit h the remanent 
magnetic field measurements in chondrules by jFu et aT] ( |2014| l 
in that the minimum midplane field strength for a wind driven 
disk at the presumed chondrule forming location (the mid¬ 
plane at 2 to 3 au) we find is smaller than the upper bound 
on the ambient field strength in the chondrule cooling reg ion 
found in tha t work. This allows the chondrules studied by jFuj 
jet al.| ( |20T4| ) to have formed in the midplane of a wind driven 
laminar disk, while also allowing for the chondrule formation 
process to locally concentrate the magnetic field. 

The influence of the traditional dead-zone picture for proto¬ 
planetary disks during various stages in pla net formation has 
been examined in some detail. For example, [Gressel, Nelson,| 
[& Turner! ( [201 1| |2012| l showed that the stirring of planetesi- 
mals by the propagation of waves from the active region into 
the dead zone can stir up small planetesimals (< 10 km) and 
cause them to undergo collisional destruction, or prevent their 
runaway growth. In recent work (Nelson et al. 2014, in prep.), 
it has furthermore been shown that low mass planets orbiting 
in disk models with traditional dead zones are likely to un¬ 
dergo very rapid inward migration because the weak stresses 
there are unable to prevent saturation of the corotation torque. 
In principle, the new emerging picture of PPDs - that is, as¬ 
suming that Maxwell stresses operate near the midplane be¬ 
cause of the Hall effect - can ameliorate these problems by 
removing the turbulent stirring, and by providing a significant 
stress that can prevent corotation torque saturation. Full sim¬ 
ulations are clearly required to examine whether or not the 
latter e ffect can be realized during full nonlinear evolution. 
Finally, [Gressel et al.[([2013|l examined gap formation and gas 
accretion onto a giant planet within a disk with a traditional 
dead zone, and observed that the gas in the gap region became 
enlivened because of enhanced penetration of X-rays and cos¬ 
mic rays. This created a turbulent and chaotic gas flow into the 
planet Hill sphere, and the enhanced magnetic coupling also 
led to the launching of a magnetocentrifugally accelerated jet 


from the circumplanetary disk. It remains unclear whether 
or not these effects will arise when ambipolar diffusion is in¬ 
cluded, given the low density of the gap region. 
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